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I. 



INTRODUCTION 



A . BACKGROUND 

Marine vehicles that require high maneuverability, such as 
hydrofoils, require quick response, control and guidance to 
maintain accurate track keeping characteristics. This is 
accomplished through successful path planning, navigation, 
guidance and autopilot design [Ref. 1]. 

Sufficient information is obtained from charted obstacles 
and the environment for smooth paths to be generated for the 
vehicle to follow [Ref. 2] . A certain level of feedback is 
provided through the use of sonar and acoustics in order to 
replan a path when uncharted obstacles are encountered or when 
the mission objectives are changed. Based on inertial 
positional information, the guidance law provides the 
appropriate vehicle heading by suitable use of the vehicle 
effectors such as rudders, dive planes, and cross body 
thrusters. However, lags in obtaining and processing inertial 
positional information can limit vehicle reaction time. In 
addition, the guidance law must be as fast as possible in 
order to ensure accurate path keeping characteristics [Ref. 
3] . Therefore, the stability of the combined guidance/control 
scheme becomes an issue that needs to be analyzed. 
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B. OBJECTIVE OF THIS THESIS 



This thesis investigates the effects of positional 
information time lags on guidance and control. The 
relationship of the critical visibility (minimum vehicle 
lookahead distance) versus the controller time constant and 
its effect on the stability of the combined guidance/control 
scheme is analyzed. Results are presented based on 
computations using time domain and frequency response 
techniques. All computations are performed for a dynamic 
model of the NPS AUV II [Refs. 4 and 5] for which a complete 
set of hydrodynamic coefficients and geometric properties is 
available . 

C. THESIS OUTLINE 

Chapter II presents the vehicle equations of motion in the 
horizontal plane, the equations that govern steering control, 
and the guidance scheme used in this research. 

Chapter III presents the stability analysis based on 
eigenvalue computation of the combined guidance and control 
scheme developed in Chapter II. Results are presented based 
on the stability curves developed for first, second and third 
order approximations for time lag in the commanded vehicle 
heading in the control law. 

Chapter IV presents an exact computation of the stability 
curves based on frequency response methods, which utilizes the 
Nyquist criterion for stability. 
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II. EQUATIONS OF MOTION, CONTROLLER DESIGN AND GUIDANCE SCHEME 



A. INTRODUCTION 

The equations of motion that describe maneuvering of a 
marine vehicle in the horizontal plane are presented in this 
chapter. A linear full state steering feedback control law is 
designed based on the three equations in sway, yaw, and rate 
of change of heading angle. The guidance scheme is based on 
pure pursuit guidance for path following along straight line 
segments . 

B. EQUATIONS OF MOTION 

Restricting our attention to motion in the horizontal 
plane (steering control), the mathematical model consists of 
the nonlinear sway and yaw equations of motion only. With a 
moving coordinate frame fixed at the vehicle's geometric 
center, Newton's equations of motion are 



m(v+ur+x G r-y G r 2 ) =Y (2.1) 



I z t+mx G {v+ur) -my G vr=N (2.2) 
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where v and r are the relative sway and yaw velocities of the 
moving vehicle with respect to the water; m is the mass of the 
vehicle; Xq, y G are the respective lateral and longitudinal 
positions of the center of gravity; and Y,N represent the 
total excitation sway force and yaw moment, respectively. 
These forces can be expressed as the sum of quadratic drag 
terms and first order memoryless polynomials in v and r, which 
represent the added mass and damping due to the vehicle's 
motion through the water. In this way, Y and N can be 
represented by 



Y = -P^y^r + ^C 3 ( Y^v+Yur) +-£{ 2 Y t 



\JLV 



_ _p 
2 



J tail y V+£r 2 



N = ■£-t 5 N f r+-P-lt i {N^v+N r ui) +-£fi 3 N v uv 
2 2 2 



where C is the vehicle length, and Y a , N a represent partial 
derivatives of Y and N with respect to a, is the drag 

coefficient, and 8 is the rudder angle. 

Equations (2.1) and (2.2) can be nondimensionalized with 
respect to the constant forward speed u, and the vehicle 
length C ; with the dimensionless time variable being tu/C . 
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The cross flow integral drag terms become important for 
hovering operations or low speed manuevering, whereas at 
relatively high speeds and low angles of attack (with respect 
to the water), the steering response is predominantly linear. 

The model becomes complete with the addition of the 
expressions for the vehicle yaw rate 



i| s-r 



(2.3) 



and the inertial position rate (horizontal plane) 



y=usim|r + vcosi|/ (2.4) 

where \\f is the heading angle as shown in Figure 1. 

Equations (2.1), (2.2), and (2.3) can be written as a set 

of three coupled linear differential equations of the form 

ijr=r (2.5a) 



v’=a 11 uv+a 12 ur+jb 1 u 2 5 



(2.5b) 



r=a 12 uv+a 22 ux+b 2 u z t> 



(2.5c) 



5 



y v 




Figure 1. Vehicle Geometry and Definition of Symbols. 
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where the coefficients a i;j and b t are functions of the vehicle 
hydrodynamic and geometric properties, \| / is the heading angle 
and 5 is the rudder angle. 

The linearized form of equation (2.4) at a nominal u is 

y=ui|f + v (2.6) 



Hence, the final set of linearized equations of motion used 
for steering control are 



i|/=r 



(2.7a) 



v=a 1:L u v+a 12 ur +1^ 



(2.7b) 



r=a 12 uv+a 22 ur+jb 2 u 2 S 



(2.7c) 



y=u\ \f + v 



(2 . 7d) 



at any nominal u. 

This system of equations becomes the basis for the 
controller design. 



7 



C. CONTROLLER DESIGN 



Equations (2.7a), (2.7b), and (2.7c) govern the steering 

control of the model used in this research. 

These equations can be represented in the form 



x-Ax+bb 



( 2 . 8 ) 



where the state vector equation is 

x = [\|r, v, r] T (2.9) 

Linear full state feedback is introduced in the form 

6=k 1 (ty—ty c ) +k 2 v+k 3 r (2.10) 

where vjl c is the commanded heading angle, and the gains k x , k 2 , 
k 3 are computed by pole placement such that the closed loop 
system of equations (2.8), (2.10) has the desired dynamics. 

The existence of the difference (y-\)/ c ) in the control law 
(2.9) has the effect of trying to point the vehicle's 
longitudinal axis towards the desired heading. 

If the desired characteristic equation is selected so that 
its time constant is t c dimensionless seconds, its general 
form becomes 
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where 



«! = 





The controller gains are computed from 



Jc x = ^ (2.12a) 

(b^-b^) u 3 



( ^ia 22 ~b 2 a 12 )u 3 k 2 +(b 2 a 11 -b ± a 21 ) u 2 k 3 = a z +b 2 u z k x (2.12b) 

b Y u 2 k 2 ^b 2 u z k^-a x - [a^a 22 ) u (2.12c) 

Due to the explicit dependence on u, the gains k 1( k 2 , k 3 are 
continuously updated for different nominal forward speeds. 

D . GUIDANCE SCHEME 

The guidance scheme used in this research is an 
orientation based command scheme. In this scheme, the 



commanded vehicle heading angle vj/ c is a function of the line 
of sight angle a between the actual vehicle position and a 
reference position on the nominal path which is located at a 
constant lookahead or visibility distance d ahead of the 
vehicle as shown in Figure 1 . 

This line of sight angle is defined by 



The autopilot is then called upon to deliver the commanded 
heading angle \j/ c . The simplest orientation guidance law is 
pure pursuit guidance where the commanded heading angle y c in 
the control law (2.9) equals the line of sight angle a in 
(2.13). 

Then \|/ c is defined by 




d 



(2.13) 



iJj c = -arctan-=f 



y 

d 



(2.14) 



For relatively small angles 



tc = -arctan-g * -■£ 



(2.15) 
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This results in a three state heading autopilot design 
which can exhibit very robust characteristics . However, the 
stability of this scheme may become an issue when transiting 
along straight line segments since the commanded heading angle 
is a function of the vehicle response, i.e., autopilot 
response is limited by time lags in vehicle dynamics. We will 
next examine the stability of this scheme. 
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III. EIGENVALUE ANALYSIS 



A. INTRODUCTION 

In this chapter, we present the stability analysis, based 
on eigenvalue computation, of the combined guidance and 
control scheme developed in Chapter II. A brief presentation 
of the stability curve is developed first for the case of zero 
time lag. Incorporation of a nonzero time lag in the 
positional information formally results in an infinite state 
space system. This is truncated into first, second, and third 
order approximations. The resulting eigenvalue problems are 
then analyzed with emphasis on stability curves. 

B. STABILITY CONSIDERATIONS AND TIME LAG 

1. Stability Considerations 

In pure pursuit guidance, where the commanded heading 
angle \j/ c in the control law equals the line of sight angle a, 
the trivial equilibrium state which corresponds to straight 
line motion is characterized by 

i|r = v=r=y=0 (3.1) 

Linearization of the state equations in the vicinity of (3.1) 
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produces the linear system 



x-Ax, 



where the complete state vector is 



(3.2) 



x- [i|r, v, r,y] T 



(3.3) 



Local stability properties of (3.1) are then established by 
the eigenvalues of A. The characteristic equation of (3.2) is 
a quartic of the form 



k*+Bk 2 +Ck 2 +Dk+E= 0, 



(3.4) 



where the coefficients B, C, D, E are functions of the vehicle 
properties, the lookahead distance d, and the controller gains 
kj , k 2 , k 3 . A pair of complex conjugate roots of (3.4) crosses 
the imaginary axis when the following conditions are met 
[Ref. 6] . 



BCD-B 2 E-D 2 =0 , 



(3.5a) 




(3.5b) 
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These two conditions (3.5a) and (3.5b) translate to 



a 1 d 2 +a 2 d+a 3 = 0 , 



(3.6a) 



d > 



■^i a 22 -^2 a i2 ^2 

b 2 a il~ b i a 21 



(3 .6b) 



where 



a i =a i a 2 -a 3 



(3.7a) 



_ (g 1 g 2 -2g 3 ) (b^-J^a^-b^ 2 

2 b 2 a il~ b i a 21 ( b 2 a il~ b i a 2l) U 1 (3 ' 7b) 

a _ -(^a 22 -b 2 a 12 -b 2 ) [b 1 g 1 +(b 1 a 22 -b 2 a 12 -b 2 )u]g 3 {3 7c) 

( b 2 a il~ b i a 21 ) 2U 

The conditions (3.5a) and (3.5b) determine the critical 
visibility d crit for stability. For d > d crit/ the equilibrium 
state (3.1) is stable which means that the control law will 
drive and keep the vehicle onto the straight line path. For 
d < d cric/ the equilibrium state loses its stability and the 
vehicle response becomes oscillatory as a result of the pair 
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of complex conjugate eigenvalues with positive real parts. 

Results of the critical visibility versus the 
controller time constant t c for zero time lag are presented in 
Figure 2. As expected, the higher values of t c (i.e., softer 
controller) require higher lookahead distances d for path 
accuracy. This agrees with existing guidelines that the 
guidance law must be sufficiently slower than the controller 
in order for the dynamics of the two not to interfere with 
each other and, therefore, guarantee stability of the combined 
scheme. Very high values of d correspond to a very slow 
guidance law with a loss in speed of response and path 
accuracy [Ref. 1]. 

2 . Time Lag 

In ocean vehicles, all states needed in the control 
law are readily available at the required rate, with the 
possible exception of the positional information y (Figure 1) . 
The latter may require a significant data analysis and 
reduction of sonar returns and inertial navigational 
information. As a result, it is likely that a time lag will 
exist in the positional information y which is used in the 
guidance law. 

With the introduction of a time lag T (sec) , the commanded 
rudder angle y c in the pursuit guidance control law becomes 
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d 



Figure 2. Critical Visibility Distance d versus t c for 
Zero Time Lag T. 
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(3.8) 






-arctan 



y(t-r) 

d 



and for relatively small angles 



Vc 



yit-T) 

d 



The resultant rudder control law (2.9) becomes 



6 +k x XiLIl + k 2 v+k 3 r 



After some algebra, 
motion (2.7a), (2.7b) 



the resultant linearized equa 
, (2.7c), (2.7d) become 






k 

v=j b 1 u 2 k 1 \lr+ (a^u+fyu 2 ^) v+ {a^u+b^u 2 ^) r+b 1 u z -^-y{ t-T) 



k 

i=b 2 u z k 1 ^/ + ( a 21 u+b 2 u z k 2 ) v+ ( a 22 u+b 2 u z k 3 ) r+b 2 u z -^-y{ t-T) 



(3.9) 



(3.10) 



ions of 



(3.11a) 



(3.11b) 



(3.11c) 
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y-u\ \r + v 



(3 .lid) 



The Taylor expansion for the term y(t-T) is 



y{t-T) » y-Ty+^ly-^-y+. . . 

2 6 



(3.12) 



We can now examine the stability of the previously developed 
control scheme for first, second and third order 
approximations for time lag in the control law. 



C. FIRST ORDER APPROXIMATION FOR TIME LAG 

For a first order approximation for time lag T 



y ( t-T) = y-Ty 



(3.13) 



where 



y=ui|r + v 



(3.14) 



After some algebra, the resultant linearized equations of 
motion (3.11a), (3.11b), (3.11c), (3. lid) become 
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(3.15a) 



i|f=r 



k £ 

v=(h 1 u 2 ^ 1 -i 5 1 u 3 -^D }f f *(a 11 u+b l u 2 k 2 -b 1 u 2 -^T) v+ (a^u+i^u 2 ^) r+i) 1 u 2 -^y (3 # 15b) 



k 

r= (i» 2 u 2 ic 1 -i) 2 u 3 -^r) i|; + (a 2 iU+l) 2 u 



2 ^ 2 _i, 2 



, 2*1 



D v+ (a 22 u+±) 2 u 2 ic 3 ) r+b 2 u 2 



k i 

~d y ( 3.15c) 



y=in|r + v 



( 3 . 15d) 



The state space form of the equations of motion is 



x-Ax 



(3.16) 



In matrix form, (3.15a), (3.15b), 



(3.15c), (3.15d) become 



* 




0 


0 


1 


0 


V 




(b^k^&u'-^T) 


k 

(a 11 u+b 1 u 2 k 2 -b 1 u 2 -^-T) 


(a^u+jbjid/Cj) 


Cb.u 2 ^) 


i 




(i 2 u 2 * 1 -i) 2 u 5 -^D 


(a 21 u*b 2 u z k 2 -b 2 u 2 ^T) 


(a 22 u+b 2 u 2 k 3 ) 


(l> 2 u 2 ^) 


y 




u 


1 


0 


0 



* 



V 

r 



(3.17) 



y 



The standard eigenvalue problem (3.17) is solved to 
determine numerically the critical d versus t c curve for a 
given value of T. Appendix A presents the program used for 
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the first order approximation for time lag. 

Figure 3 shows the resulting stability curves for the 
first order case. It shows the relation between the critical 
visibility d and the controller time constant t c for values of 
time lag of 0.0, 0.5, 1.0, 1.5 and 2.0 seconds. 

As expected, the stability curve shifts to the right with 
increasing values of time lag. For the same time constant t c , 
the critical visibility distance increases with corresponding 
increases in the amount of time lag. 

D. SECOND ORDER APPROXIMATION FOR TIME LAG 

For a second order approximation for time lag T, we have 



y ( t-T) = y-Ty+^-y 



(3.18) 



where 



y=m|r + v, y-ur+v (3.14), (3.19) 

Following some algebra, the resultant linearized equations of 
motion (3.11a), (3.11b), (3.11c), (3. lid) become 
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d 

Figure 3. First Order Approximation For Time Lag: 

Critical Visibility Distance d versus t c for 
Time Lags T = 0.0, 0.5, 1.0, 1.5, 2.0 sec. 
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(3.20a) 



where 



^=r 



+A 22 v+A 22 r+A 24 y (3.20b) 



r=A 31 i|r+.A 32 v+.A 33 r+.A 34 y (3.20c) 



y=ux|/ + v (3 . 20d) 



A 2 i 



-jb 3 u 3 — t T 



__1 

d 



l-b r u 2 -^-T 2 
1 2d 



■^22 



a n u+ -^i u 2 k 2 ~b 1 u 2 T 

l-bu 2 -^ 2 
1 2d 



2 3 



a i2 U+b l u2k 2 +b l u3 ~2^ T2 

1-b u 2 ^-T 2 
1 2d 



22 



2*1 



■*24 * 



w-i 

l-b.u^T* 



k k {b 1 u 2 k 1 -b 2 u 2 -^rT) 

^ 3 i = {b 2 u z k x -b 2 u 2 -^T) ^{b 2 u 2 ^T 2 ) 



32 



= (a 21 u+b,u z kj-byU 2 — ^ T ] + (b 2 u 2 -^T 2 ) 



k i 



k i 



k 

( a n u +b 1 u z k 2 -b r u z —±T) 



2 d 



(1 -b x u 2 ^T 2 ) 



^33 = 



(a 22 u+b 2 u 2k 3 +b z u 3 -^±T 2 ) + (b 2 u 2 ^bT 2 ) 



(a 
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u +b 1 u z k 2 +jb t u 3 r 2 ) 



(l-^u 2 — ^T 2 ) 
1 2d 



^34 = 



(b 2 u 2 ^) + (b 2 u 2 -^T 2 ) 
2 d 2 2d 



, 2 k i 






(i -d 1 u 2 -^i'r 2 ) 



As with the first order case, the state space equations of 
motion reduce to 
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x-Ax 



(3.16) 



In matrix form, (3.20a), (3.20b), (3.20c), (3 . 2 Od ) 







0 0 10 




t 


V 




■^21 -^22 -^23 -^24 




V 


r 




-^31 As2 A33 ^34 




I 


y 




u 1 0 0 




y 



become 



(3.21) 



Simular to the first order case, the standard eigenvalue 
problem (3.21) is solved using Fortran programming and MATLAB 
in order to develop the stability characteristics of the 
second order case. Appendix B presents the program used to 
develop the stability characteristics for this case. 

Figure 4 shows the resulting stability curves for this 
case. The results of the critical visibility d versus the 
controller time constant t c for time lags of 0.0, 0.5, 1.0, 
1.5 and 2.0 secs are shown. 

As with the first order case, the higher values of t c 
require higher lookahead distances d for path accuracy. High 
values of d correspond to a slow guidance law with a loss in 
speed of response and path accuracy, and the stability curve 
shifts to the right with increasing values of time lag. 
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d 



Figure 4. Second Order Approximation For Time Lag: 

Critical Visibility Distance d versus t c for 
Time Lags T = 0.0, 0.5, 1.0, 1.5, 2.0 sec. 
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However, due to being more accurate than the first order case, 
the second order approximation provides tighter stability 
curves for the same time lags. 



E. THIRD ORDER APPROXIMATION FOR TIME LAG 

For a third order approximation for time lag T 



yU~T) = y-Ty+I-y-I-y 



(3.22) 



where 



y=ux J/ + v, y=ur + v, y=ui + v (3.14), (3.19), (3.23) 

After a considerable amount of algebra and following a simular 
procedure as with both the previous two cases, the resultant 
linearized equations of motion reduce to the state space form 

Bx=Ax (3.24) 



or 



x=B 1 Ax 
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The resultant linearized equations of motion (3.11a), 



(3.11b), (3.11c), (3. lid) become 



i|;=r 



B 33 i+B 3S^2 = A 31$ +A 32 V l +A 33 V 2 +A 34 r+A 3S V 



B 53^ +B 55^2 = A ^ +A S2 V l +A 53 V 2 +A 5i r+A S5V 



y=u^ + v 1 



= 



v. 



where 



B 



33 



= jb.U 3 — ^T 3 

1 6 d 



B 



35 




A 



31 



J b 1 u 2 A- 1 -i) 1 u 3 



d 



T 



(3.25a) 



(3.25b) 



(3.25c) 



(3 . 25d) 



(3 . 25e) 
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A 32 = a^u+b^kz-b^u 2 -±T 



a z3 = ai2 u+b i u2k i +b i u3 



A 3< 




3 5 



= jb,u 2 -^r 2 -i 

1 2d 



and 



’ ltb > ul Td T> 



B 



55 




*i. 



A>1 = b 2 uZ - b 2 u3 -J T 



A 52 = a 21 u+b 2 u 2 k 2 -b 2 u 2 -±T 
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^53 = a 22 U+b 2 U 2 k 3 + b 2 U 3 T 2 



^54 = b 2 u2 ^J 



A„ = b y u 2 ^-T 2 
55 2 2d 



In matrix form, (3.25a), (3.25b), (3.25c), (3.25d), (3.25e) 

become 



1 


0 


0 


0 


0 




* 




0 


0 


1 


0 


0 




* 


0 


1 


0 


0 


0 








0 


0 


0 


0 


1 




v i 


0 


0 


*33 


0 


*35 




r 


= 


*31 


*32 


*33 


*34 


*35 




r 


0 


0 


0 


1 


0 




y 




u 


1 


0 


0 


0 




y 


0 


0 


*53 


0 


*55 




*2 




*51 


*52 


*53 


*54 


*55 




V 2 



(3.26) 



Simular to both the first order and second order cases, 
the generalized eigenvalue problem (3.21) is solved using 
Fortran programming and MATLAB in order to develop the 
stability characteristics for the third order case. Appendix 
C presents the program used to develop the stability 
characteristics for this case. 

Figure 5 shows the resulting stability curves for this 
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Figure 5. Third Order Approximation For Time Lag: 

Critical Visibility Distance d versus t c for 
Time Lags T = 0.0, 0.5, 1.0, 1.5, 2.0 sec. 
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case. The results of the critical visibility d versus the 
controller time constant t c for time lags of 0.0, 0.5, 1.0, 
1.5 and 2.0 seconds are shown. As with the first and second 
order cases, the higher values of t c require higher lookahead 
distances d for path accuracy, high values of d correspond to 
a slow guidance law with a loss in speed of response and path 
accuracy, and the stability curve shifts to the right with 
increasing values of time lag. However, - due to being more 
accurate than either of the other two cases, the third order 
approximation provides the tightest stability curves for the 
same time lags. 



F. COMPARISON OF RESULTS 

Results of the critical visibility d versus the controller 
time constant t c for first, second and third order 
approximations for time lag have been obtained. Figures 3, 4 
and 5 presented the stability curves for first, second and 
third order approximations for time lags of 0.0, 0.5, 1.0, 1.5 
and 2.0 secs. These figures have shown that stability of the 
control scheme decreases with increasing time lag, and for the 
same time constant, the critical visibility distance increases 
with increasing time lag. 

Figures 6, 7, 8 and 9 present a comparison of first, 
second and third order approximations for each time lag 
considered. These figures show that for each time lag, the 
third order model presents the largest region of stability for 
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d 

Figure 6. Critical Visibility Distance d versus t c for 
First, Second and Third Order Approximation 
Cases and Time Lag T = 0.5 sec . 
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Figure 7. Critical Visibility Distance d versus t c for 
First, Second and Third Order Approximation 
Cases and Time Lag T = 1.0 sec . 
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Figure 8. Critical Visibility Distance d versus t c for 
First, Second and Third Order Approximation 
Cases and Time Lag T = 1.5 sec. 
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d 

Figure 9. Critical Visibility Distance d versus t c for 
First, Second and Third Order Approximation 
Cases and Time Lag T = 2.0 sec. 
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vehicle straight line path accuracy. In addition, for the 
same time constant, the critical visibility is least for the 
third order case. This demonstrates that the minimum 
lookahead distance required for straight line motion stability 
decreases with controller time lag accuracy. 

It can be seen that the first order approximation is the 
most conservative for design purposes since it predicts the 
highest required minimum visibility distance. As expected, 
the differences among the three approximations are more 
pronounced as the time lag increases and as the controller 
time constant becomes smaller; i.e., tighter control. 

Figures 10, 11, 12 and 13 present the differences in 
critical visibilities among first, second and third order 
approximations for each individual time lag considered. A 
comparison of these figures demonstrates that the difference 
in critical visibility between respective orders (time lag 
models) decreases with decreasing time lag. This indicates 
that for low time lags, controller time lag accuracy is not as 
crucial for maintaining vehicle straight line path accuracy as 
that of higher time lags. The greater the time lag, the 
greater the difference in critical visibilities among time lag 
models, and the more crucial is the controller time lag 
accuracy for stability. 

In general, if the guidance law is slow enough to allow 
sufficient time for the controller to react, stability of 
straight line motion is guaranteed [Ref. 1]. 
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dT-dO 

Figure 10. Difference of Critical Visibility Distances at 
Time Lags of T = 0.5 and 0 . 0 sec versus t c for 
First, Second and Third Order Cases. 
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dT-dO 

Figure 11. Difference of Critical Visibility Distances at 
Time Lags of T = 1.0 and 0.0 sec versus t c for 
First, Second and Third Order Cases. 
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dT-dO 



Figure 12. Difference of Critical Visibility Distances at 
Time Lags of T = 1.5 and 0 . 0 sec versus t c for 
First, Second and Third Order Cases. 
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dT-dO 

Figure 13. Difference of Critical Visibility Distances at 
Time Lags of T = 2.0 and 0.0 sec versus t c for 
First, Second and Third Order Cases. 
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G. NORMALIZATION OF RESULTS 



In an attempt to characterize the stability analysis 
results for all values of time lag, we proceed as follows: If 
the vehicle side slip velocity v is very small so that it can 
be neglected, the equations of motion (2.7a), (2.7b), (2.7c), 
and (2.7d), take the form 



\|/=r 



r=a 22 ur+b 2 u 2 Z> 



y=ui\i 



since v = 0 . The control law is 

b=k 1 (i|f-t|/ c ) +k 3 r, 

where the first order approximation for the commanded heading 
angle is 



*c = 



yzly 

d ' 



Based on the above equations, we can get the characteristic 
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equation of the system as 



. 3 _ 



(a 22 u+b 2 k 2 u 2 ) s 2 - (b 2 k 1 u 2 -b 2 k 2 u 3 T^) s-b 2 k 2 u 2 = 0. 



Applying Routh's criterion to this cubic equation, we find 
that loss of stability occurs when 



{a 22 u+b 2 k 2 u 2 ) {b 2 k x u 2 -b 2 k r u 3 T-i ) 



~b 2 k iu 3 



1 

d' 



from which we can get the critical visibility distance 



d 



1 

a 22 +b 2 k Z U 



+ Tu 



Equation (3.27) can be written as 



(3.27) 



b T d 0 

T 



u, 



(3.28) 



where d 0 is the critical visibility distance for T = 0, and d T 
the critical visibility distance for a nonzero T. In 
dimensionless quantities, equation (3.28) is 
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dr-dp _ 



1 



(3.29) 



T 



where u = 2.0 ft/sec is the vehicle speed, and t = 7.3 ft is 
the vehicle length. 

A comparison among the approximate expression (3.29) and 
the first, second, and third order approximations is shown in 
figures 14, 15 and 16. The agreement is better for the first 
order approximation, as expected, and also for the higher 
controller time constant t c . This is because a high t c results 
in soft vehicle response with negligible amounts of side slip. 
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(dT-dO)/T 



Figure 14. Normalized Critical Visibility versus t c for 
First Order Approximation For Time Lags 
T = 0.5, 1.0, 1.5, 2.0 secs. 
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Figure 15 . Normalized Critical Visibility versus t c for 
Second Order Approximation For Time Lags 
T = 0 . 5 , 1 . 0 , 1 . 5 , 2.0 secs. 
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Figure 16. Normalized Critical Visibility versus t c for 
Third Order Approximation For Time Lags 
T = 0.5, 1.0, 1.5, 2.0 secs. 
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IV. FREQUENCY RESPONSE ANALYSIS 



A. INTRODUCTION 

The previous analysis using Taylor expansions of y(t-T) 
breaks down for approximations beyond third order. The reason 
for this is that for higher than third order expansions, the 
matrix B in the generalized eigenvalue problem (3.23) becomes 
singular. Therefore, it was felt that a different technique 
should be sought in order to obtain an exact computation of 
the stability curves and also to check the validity of our 
calculations. In this chapter, we present a technique based 
on frequency response methods, which utilizes the Nyquist 
criterion for stability. The resulting equation is solved 
using Newton's iteration. 

B . NYQUIST CRITERION 

From our previously developed model, the linearized 
horizontal equations of motion are 



i|r=r 



(2.7a) 



v=a 11 uv+a 12 ur+jb 1 u 2 6 



(2.7b) 



47 



(2.7c) 



r=a 12 uv+a 2Z uz+b 2 u 2 !> 



y-u\ ]f + v 



(2. Id) 



at any nominal u. The control law is 



t>=k x (4f-ilf c ) +k 2 v+k 2 r 



( 2 . 10 ) 



where the commanded heading angle is 



V. 



y(t-T) 



(3.9) 



The Laplace transform of equation (3.9) is 



Yc = 




(4.1) 



and the resultant rudder control law (2.10) becomes 



6 = k x tj/ + k 2 v+ k 3 r + 222. y e ~ T3 



d 



(4.2) 



Following some algebra, the resultant linearized equations of 
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motion (2.7a), (2.7b), (2.7c), (2.7d) become 



i| ;=r (4 . 3a ) 



v= (jb 1 u 2 ic 1 -±) 1 u 3 -^'r) i|»+ (a 11 u+h 1 u 2 ic 2 -jb 1 u 2 ^r) v+ (a^u+^u 2 ^) r+b 1 u 2 -^ye' T3 (4.3b) 



r= (i) 2 u 2 ic 1 -£i 2 u s -^r) i( >+ (a 21 u+b 2 u 2 k 2 -b 2 u 2 -^T) v+ (a 22 u+b 2 u 2 k 3 ) r*b 2 u 2 -^ye' Ts (4 . 3c ) 



y=uijr + v 



(4.3d) 



Equations (4.3a), (4.3b), (4.3c), (4.3d) reduce to the state 

space form 



x=Ax 



(3.16) 



The matrix form becomes: 



♦ 




0 


0 


1 


0 






(^u 2 ^) 


(a 11 u+ J b 1 u 2 /c 2 ) 


[a 12 u+b 2 u 2 k 2 ) 


( J b 1 u 2 -^e' rs ) 


r 




{b 2 u 2 k 2 ) 


(a 21 u+b 2 u z k 2 ) 


(a 2 2 u +b 2 u 2 k 2 ) 


(i5 2 u 2 -^e' rs ) 


y 




u 


1 


0 


0 



* 



v 

r 



(4.4) 



y 
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The characteristic equation of the form [A-Is] = 0 becomes 
after a considerable amount of algebra 

s*+A 3 s 2 +A 2 s 2 +A 1 s+ {B 2 s 2 +B 1 s+B 0 ) De~ Ts = 0 (4.5) 

where 

^3 = -( a n +a 2 2 ) u -{b^+bzkj u 2 

^2 ~ ^ a ii a 22 -a i2 a 21^ UZ+ ^ a il^2~ a 21^- > l^ U + ^ a 22^" ) l _a i2^ ) 2 ^ U ^ ^2 ~^2 U 2 

■^i - ~ ( a i2^2 _a 22^i ) ^ k-JDG — ( 3 2 i-bi ~3 11 i ) 2 ) ^ k^ 



B 2 = ~b 1 u 2 k 1 



B 1 = b 2 u i k 1 



B 0 = (a u i 2 -a 2lJ b 1 ) u 4 Jc, 



and 



D = 



1 

d 
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The characteristic equation of the system is written as 



1 + K G ( s) = 0 



(4.6) 



where 



G{s) = (B 2 s 2 +B 1 s+B 0 ) e~ Ts 

s(s 2 +A 1 s z +A 2 s+A 1 ) 



(4.7) 



is the transfer function, and we have denoted 

K = D 



(4.8) 



With s=jco, the phase angle $ is represented by 



<j> = ZG(jg>) 



(4.9) 



and it follows that 



<t> = Z(e" TJ “) - Z(jg>) + Z ( -S 2 o) 2 +B 1 jg)+S 0 ) - Z (-j(o 3 -A 3 o) 2 +A 2 jco+A 1 ) 



<J ) = -uT - - + tan 1 ( 



B 1 0) 



■B 0 -S 2 6) 2 



) - tan' 1 ( (4.10) 

- 2 1 3 0) 2 +A 1 
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For stability, the Nyquist criterion states (Figure 17) that 
at the solution of the equation 



<J> ( a> ) = -7i 



(4.11) 



which is phase cross over frequency co = o^, the gain margin 
must be equal to 1 



|iCG(jco 1 ) | = 1 



(4.12) 



where K = D = 1/d. It follows that equation (4.11) becomes 



-to, T-— +tan~ 1 ( 
1 2 



B i<*i 

B 0 -B 2 u\ 



-tan' 1 ( 



-( 0 ^^ 2 ( 0 1 

-A 3 (j>1+A 1 



-7i (4.13) 



Rearranging, we get 



tan 1 ( 



B i<*i 

Bq-B 2 of 



-tan 1 ( 



-AjCOi+Al 



-71 +— +0), T 

2 1 



Taking the tangent of both sides and using the identity 



tan [x-y) 



_ tan (x) -tan (y) 
1+tan (x) tan (y) 
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IV 




Figure 17. Example of the Nyquist Stability Criterion For 
Three Different Values of Gain [Ref. 7]. 
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and rearranging following some algebra, the result becomes 



B 1 o > 1 (A 1 -A 3 (Oi) - (B 0 -B 2 (Oi) (-(Oi+A 2 (i) 1 ) __ i 

(B 0 -B 2 <j>1) (- 4 , 0 ?+^) +fl 1 u 1 (-a)i+A 2 o) 1 ) tan (o) 1 T) 



From equation (4.12), it follows that 



1 

d io 1 1 -j(jil-A 2 a>l+A 2 ju) 1 +A 1 1 



Solving for d, we get 



v /(B 0 -g 2 0)^) 2 +B 1 2 (0? 
o> lV / (^-^0)1) 2 + (A,o> 3 -g>i) 2 



With a time lag T = 0, we get from equation (4.14) 



(B 0 -B 2 <j) 1) (AjGJi+j^) + 5 X 0)3 ( -Ui+.A^) = 0 



Rearranging, we get 



(B^-B^ w* + (B^-B^-BqAj) 0)3 + 5 0 ^3 = 0 



(4.14) 



(4.15) 



(4.16) 
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Equation (4.16) can be solved exactly for 0^, and then d 
can be computed from (4.15). This zero time lag calculation 



was performed in order to validate the results 
for the zero time lag solution of equation (3.6) 
using the two different techniques were found to 



C . ALGORITHM DESCRIPTION 

Introducing the terms 



«i = b 2 a z -b. 



-® 1^2 -® 2^1 -® 0^3 



a 3 -® 0^1 



and 



Pi = 



Pi - sa-ba 



P 2 B 0 + B z A 2 B x A 2 



of d versus t c 
The results 
be identical . 
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Then equation (4.12) becomes 



P 1 (0 5 + P 2 G) 3 +P 3 0) _ 1 

a 1 (i) 4 +a 2 o) 2 +a 3 tan(a>D 



(4.17) 



Rearranging, we get 



(P 1 co 5 + P 2 (i) 3 + P 3 (i)) tan (or) + a 1 co 4 +a 2 co 2 + a 3 = 0 (4.18) 



Equation (4.18) is now in the form 



f (o) = 0 



(4.19) 



Following Newton's iteration for co, we have 






j (<■>*-!> 

f '^k- 1 ) 



(4.20) 



From equation (4.20), 



f( o) = ( p 1 (i> 5 +p 2 (i> 3 +p 3 G) ) sin(or) 
+ (a 1 o> 4 +a 2 Ci) 2 +a 3 ) cos (or) 



(4.21) 
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and 



= (5P 1 ( J ) 4 +3P 2 o) 2 + P 3 ) sin (or) 

+ (p, (i) 5 +p 7 (o 2 + P-,ci)) Tcos(u>T) 

(4.22) 

+ (4a 1 cjj 3 +2a 2 o)) cos (car) 

- (a 1 o> 4 +a 2 (i) 2 +a 3 ) Tsin (coD 

Fortran programming and MATLAB were used in order to 
develop the stability characteristics for the Newton's 
iteration method. Appendix D presents the program used to 
develop the stability characteristics for this case. 

Figure 18 shows the resulting stability curves for the 
Newton's iteration method. The results of the critical 
visibility d versus the controller time constant t c for time 
lags of 0.0, 0.5, 1.0, 1.5, and 2.0 secs are shown. 

As with the first, second and third order approximation 
cases, the higher values of t c require higher lookahead 
distances d for path accuracy, high values of d correspond to 
a slow guidance law with loss in speed of response and path 
accuracy, and the stability curves shift to the right with 
increasing values of time lag. However, since the Newton's 
iteration method presents the exact solution of the frequency 
response to equation (4.19), the resulting stability curves 
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d 

Figure 18. Newton's Iteration Method: Critical Visibility 
Distance d versus t c for Time Lags T = 0.0, 
0.5, 1.0, 2.0 secs . 
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represent the precise stability characteristics for these time 
lags . 

D. RESULTS AND DISCUSSION 

Results of the critical visibility d versus the controller 
time constant t c for the Newton's iteration method has been 
demonstrated. Figure 18 presented the precise stability 
curves for time lags of 0.0, 0.5, 1.0, 1.5 and 2.0 secs. 

This figure, as well as those for the first, second and 
third order approximation cases, demonstrates that the 
stability of the control scheme decreases with increase in 
time lag and for the same time constant, critical visibility 
increases with increase in time lag. 

Figures 19, 20, 21, 22 and 23 present a comparison of 
generated stability curves among the Newton's iteration method 
and those for the first, second and third order approximation 
cases for time lags of 0.0, 0.5, 1.0, 1.5 and 2.0 secs. It 
can be seen that there is barely a slight difference in the 
stability curves for the Newton's iteration method case and 
the third order approximation case. This indicates that a 
third order approximation for time lag presents the best model 
for accounting for a positional information time lag in the 
vehicle control law. 
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d 

Figure 19. Critical Visibility Distance d versus t c for 
Newton's Iteration Method and First, Second and 
Third Order Approximations; T = 0.0 sec. 
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d 



Figure 20. Critical Visibility Distance d versus t c for 
Newton's Iteration Method and First, Second and 
Third Order Approximations; T = 0.5 sec. 
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d 



Figure 21 . Critical Visibility Distance d versus t c for 
Newton's Iteration Method and First, Second and 
Third Order Approximations; T = 1.0 sec. 
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Figure 22. Critical Visibility Distance d versus t c for 
Newton's Iteration Method and First, Second and 
Third Order Approximations; T = 1.5 sec. 
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d vs tc for 1st, 2nd, 3rd Ord & newtons method time lag T=2.0 




d 



Figure 23. Critical Visibility Distance d versus t c for 
Newton's Iteration Method and First, Second and 
Third Order Approximations; T = 2.0 sec. 
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CONCLUSIONS AND RECOMMENDATIONS 



A. CONCLUSIONS 

A methodology for considering positional information time 
lags in the control law for vehicle maneuvering in the 
horizontal plane has been presented. The relationship of the 
critical visibility versus the controller time constant for 
both zero and non-zero time lags has been established. 

Time lags have been approximated by first, second and 
third order Taylor expansions in the commanded vehicle heading 
in the control law. A comparison of the resulting stability 
curves has demonstrated that the third order approximation 
presents the greatest stability and the least critical 
visibility for a given time lag. It was found that the first 
order approximation was the most conservative for controller 
design purposes since it predicted the highest required 
critical visibility distance. In addition, the differences 
among the three approximations became more pronounced with 
increasing time lag and decreasing controller time constant 
(tighter control) . 

Further analysis was conducted to characterize the 
stability results for any value of time lag using a 
normalization procedure and the first order approximation for 
time lag in the commanded vehicle heading in the control law. 
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A comparison of the normalized stability curves for each time 
lag showed that the greatest differences occur with increasing 
order and decreasing controller time constant. This is due to 
a tight (quick) vehicle response with significant side slip at 
lower time constants . 

Frequency response techniques based on the Nyquist 
stability criterion were used to produce the exact stability 
curves for comparison with the previously generated Taylor 
expansion stability curves. These results virtually matched 
those obtained for the third order approximation time lag case 
and validated our calculations and stability curves generated 
for the case of zero time lag. These results indicated that 
the third order approximation best represents a positional 
information time lag in the control law, however for design 
purposes, the first order approximation provides the most 
conservative approach. 

B . RECOMMENDATIONS 

Some recommendations for further research are as follows : 

* Experimental verification using the NPS AUV II. 

* Consideration of time lag in the simulation of an Inertial 
Navigation System required for positional updates. 

* Analysis of positional information time lags in the motion 
stability in the vertical plane, along curved paths, or 
with respect to other guidance schemes. 
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on n o o o 



APPENDIX A 



PROGRAM THES I S 1 . FOR (Time Delay-lst Order Approx) 
BIFURCATION ANALYSIS 

IMPLICIT DOUBLE PRECISION (A-H,0-Z) 

DOUBLE PRECISION K 1 , K2 , K 3 , L , NR , NV , NDRS , NDRB , I Z , MAS S , 
£, NRDOT , NVDOT 

DIMENSION A( 4,4 ) , FV1 ( 4 ) , IV1< 4 ) ,ZZZ( 4, 4 ) ,WR( 4 ) ,WI ( 4 ) 

OPEN (11, FILE-'BIFl .RES' , STATUS= ' NEW ' ) 

OPEN ( 12 , FI LE- ' BI F2 . RES ' , STATUS = ' NEW' ) 

OPEN ( 13 , FILE-' BIF3 . RES ' , STATUS-' NEW' ) 

Vehicle Parameters: 

WEIGHT-4 35.0 



IZ 


-45.0 


L 


= 7.3 


RHO 


= 1.94 


G 


-32.2 


XG 


=0.0104 


MASS 


=WE I GHT/G 



YRDOT — 0. 00178*0. 5*RHO*L**4 
YVDOT =-0 . 034 30*0 . 5*RHO*L**3 
YR -+0. 01187*0. 5*RHO*L**3 
YV --0. 03896*0. 5*RHO*L**2 
YDRS - + 0. 02345*0. 5*RHO*L**2 
YDRB =+0. 02345*0. 5*RHO*L**2 
NRDOT — 0. 00047*0. 5*RHO*L**5 
NVDOT --0 . 00178*0 . 5*RHO*L**4 
NR --0. 01022*0. 5*RHO*L**4 
NV =-0. 00769*0. 5*RHO*L**3 
NDRS =-0.337*0. 02345*0. 5* RHO*L** 3 
NDRB -+0.283*0.02345*0. 5*RHO*L**3 
C 

WRITE (*,1001) 

READ (*,*) TCM I N , TCMAX , I TC 

WRITE (*,1002) 

READ (*,*) XDM I N , XDMAX , I XD 

XDMI N=XDMIN*L 
XDMAX=XDMAX*L 
WRITE (*,1003) 

READ (*,*) U 

C 

WRITE ( * , 1100 ) 

READ (*,*) TL 
C 

DH =( IZ-NRDOT) * ( MASS- YVDOT) - 
& ( MASS *XG- YRDOT) * ( MASS *XG-NVDOT ) 

AAl 1— ( ( IZ-NRDOT) * YV- ( MASS *XG- YRDOT ) *NV )/DH 
AA1 2= ( ( IZ-NRDOT) * ( -MASS+YR)- 
& ( MASS *XG- YRDOT ) * ( -MASS *XG + NR ) )/DH 

AA2 1= ( ( MASS- YVDOT ) *NV-( MASS *XG-NVDOT) *YV )/DH 
AA2 2 = ( ( MASS -YVDOT ) * ( -MAS S * XG+NR ) - 
& ( MASS *XG-NVDOT ) * ( -MASS + YR ) ) /DH 

BB 1 1- ( ( IZ-NRDOT) * YDRS- ( MASS* XG- YRDOT) *NDRS )/DH 
BBl 2= ( ( IZ-NRDOT) * YDRB- ( MASS* XG- YRDOT) * NDRB ) /DH 
BB2 1= ( ( MASS -YVDOT) *NDRS- ( MASS *XG-NVDOT ) *YDRS )/DH 
BB22=( (MASS-YVDOT) *NDRB-( MASS*XG- NVDOT) * YDRB ) /DH 
C 

WRITE (*,1004) 
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c 



READ (*,*) 



RATIO 



C 



C 



C 



C 



C 



C 



BB 1 =BB 1 1 + RAT 10* BB 1 2 
BB2=BB2 l+RATIO*BB22 

EPS = 1 . D-5 
I LMAX= 1500 

DO 1 1=1 , ITC 

WRITE (*,2001) I , ITC 

TOTCMIN+ ( 1-1 ) * ( TCMAX-TCMIN )/( ITC-1 ) 

TCD=TC * L/U 
POLEC=l . 0/TCD 
ADl-3 . 0 * POLEC 
AD2 = 3 . 0 * POLEC* * 2 
AD3=POLEC* * 3 
A 1 = BB 1 * U * U 
B 1 =BB2 *U *U 

Cl=-ADl- ( AA1 1 +AA2 2 ) *U 
A2 = ( BBl * AA 2 2 - B B 2 * AA 1 2 ) *U* * 3 
B 2 = ( B B 2 * AA 1 1 - B B 1 * AA 2 1 ) *U**3 
Kl=AD3/( ( BB2 * AA1 1-BB 1 * AA2 1 ) *U* *3 ) 

C2=AD2- ( AA1 1 *AA2 2-AA1 2 *AA2 1 )*U**2 + BB2*U*U‘Kl 
K2=(Cl*B2-C2*Bl )/(Al*B2-A2*Bl ) 

K3=(C2*Al-Cl*A2)/(Al*B2-A2*Bl ) 

DO 2 J=1 , IXD 

XD=XDMIN+ ( J-l ) * ( XDMAX-XDMIN )/( IXD-1 ) 

CALL L INEAR ( K1 , K2 , K3 , U , XD , AAl 1 , AA12 , AA21 ,AA22,BBl,BB2,A,TL) 
CALL RG( 4 , 4 ,A,WR,WI ,0 , ZZZ, I VI , FV1 , I ERR) 

CALL DSTABL ( DEOS ,WR,WI , FREQ) 

IF ( J.GT. 1) GO TO 10 
DEOS 00= DEOS 
XDOO =XD 
L L = 0 
GO TO 2 

10 DEOSNN-DEOS 

XDNN =XD 
PR=DEOSNN * DEOSOO 
IF (PR.GT.O.DO) GO TO 3 
LL=LL+ 1 

IF (LL.GT.3) STOP 1000 
I L = 0 

XDO=XDOO 
XDN=XDNN 
D EOS 0= DEOSOO 
DEOSN-DEOSNN 
6 XDL=XDO 

XDR=XDN 
DEOSL=DEOSO 
DEOSR=DEOSN 
XD= ( XDL+XDR ) /2 .DO 

CALL LI NEAR ( Kl,K2,K3,U,XD,AAll,AAl2,AA21,AA22,BBl,BB2,A,TL) 
CALL RG ( 4, 4,A,WR,WI , 0 , ZZ Z , IVl , FVl , I ERR) 

CALL DSTABL( DEOS ,WR,WI , FREQ) 

DEOSM=DEOS 

XDM=XD 
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PRL=DEOSL * DEOSM 

PRR=DEOSR*DEOSM 

IF ( PRL.GT. 0.D0) GO TO 5 

XDO-XDL 

XDN=XDM 

DEOSO=DEOS L 

DEOSN= DEOSM 

I L= I L+ 1 

IF ( IL.GT. ILMAX) STOP 3100 
DI F=DABS ( XDL-XDM ) 

IF (DIF.GT.EPS) GO TO 6 

XD-XDM 

GO TO 4 

5 IF ( PRR.GT. 0 .DO ) STOP 3200 

XDO=XDM 
XDN-XDR 
DEOSO-DEOSM 
DEOSN-DEOSR 
I L=I L+l 

IF ( IL.GT. ILMAX) STOP 3100 
DI F=DABS ( XDM-XDR ) 

IF (DIF.GT.EPS) GO TO 6 
XD-XDM 

4 LLL=10+LL 

WRITE ( LLL , * ) XD/L , TC 
3 XDOO-XDNN 

DEOSOO-DEOSNN 
2 CONTINUE 
1 CONTINUE 



1001 


FORMAT 


( ' 


ENTER 


MIN, 


MAX, 


AND 


INCREMENTS 


OF 


TC' ) 


1002 


FORMAT 


( ' 


ENTER 


MIN, 


MAX, 


AND 


INCREMENTS 


OF 


XD' ) 


1003 


FORMAT 


( ' 


ENTER 


U' ) 












1004 


FORMAT 


( * 


ENTER 


BOW/STERN 


RUDDER RATIO' ) 






1100 


FORMAT 


( ' 


ENTER 


TIME 


LAG 


TL' ) 








2001 


FORMAT 


(2 


15) 















END 

C 

SUBROUTINE DSTABL ( DEOS , WR , WI , OMEGA) 

IMPLICIT DOUBLE PRECISION (A-H,0-Z) 

DIMENSION WR( 4 ) ,WI ( 4 ) 

DEOS--1 . 0D+20 
DO 1 1-1,4 

IF (WR( I) .LT.DEOS) GO TO 1 
DEOS-WR ( I ) 

I J= I 

1 CONTINUE 

OMEGA-WI ( I J ) 

OMEGA-DABS ( OMEGA ) 

RETURN 

END 

C 

SUBROUTINE LI NEAR ( Kl , K2,K3,U,XD, AAl 1 , AAl 2 , AA2 1 , AA2 2,BBl,BB2,A,TL) 
IMPLICIT DOUBLE PRECISION (A-H.O-Z) 

DOUBLE PRECISION K1,K2,K3 
DIMENSION A( 4 , 4 ) 

A( 1 , 1 ) -0 . 0D0 
A ( 1 , 2 ) -0 . 0D0 
A ( 1 , 3 ) = 1 . 0D0 
A( 1 , 4 ) -0 . 0D0 
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A(2,l)= BBl *U*U*K1 * ( XD-U*TL )/XD 

A(2,2)=U*(BB1*U*( K2 *XD-K1 *TL) +AA1 1 *XD)/XD 

A ( 2 , 3 )=AAl2*U + BBl *U*U*K3 

A ( 2 , 4 } = BBl *U*U*Kl/XD 

A(3,l)= BB2*U*U*Kl*(XD-U*TL)/XD 

A(3,2)=U*(BB2*U*( K2 * XD-K1 *TL ) +AA2 1 *XD ) /XD 

A( 3, 3 )=AA22*U+BB2*U*U*K3 

A( 3,4)** BB2*U*U*Kl/XD 

A ( 4,1) =U 

A ( 4 , 2 ) = 1 . ODO 

A( 4, 3 ) “0 . ODO 

A ( 4 , 4 ) = 0 . ODO 

RETURN 

END 



70 



on o non 



APPENDIX B 



PROGRAM THES I S2 . FOR (Time Delay-2nd Otdet Approx) 
BIFURCATION ANALYSIS 

IMPLICIT DOUBLE PRECISION (A-H.O-Z) 

DOUBLE PRECISION Kl , K2 , K 3 , L , NR , NV , NDRS , NDRB , I Z , MASS , 
& NRDOT , NVDOT 

DIMENSION A{ 4,4), FVl( 4) , IVl ( 4 ) , ZZZ( 4 ,4),WR(4),WI(4) 

OPEN (ll,FILE='BIFl.RES' , STATUS = ' NEW ' ) 

OPEN ( 12 , FILE-' BI F2 . RES' , STATUS-' NEW' ) 

OPEN < 13 , FILE-' BI F3 . RES' , STATUS- 'NEW' ) 

Vehicle Parameters: 

WEIGHT-435.0 
IZ =45.0 
L =7.3 

RHO =1.94 
G =32.2 

XG =0.0104 
MASS =WEIGHT/G 

YRDOT — 0.00178*0. 5*RHO*L**4 
YVDOT =-0. 03430*0. 5*RHO*L**3 
YR =+0. 01187*0. 5*RHO*L**3 
YV =-0 . 03896*0 . 5*RHO*L**2 
YDRS -+0. 02345*0. 5*RHO*L**2 
YDRB =+0 . 02345*0 . 5*RHO*L**2 
NRDOT =-0. 00047*0. 5*RHO*L**5 
NVDOT — 0 . 00178*0 . 5*RHO*L**4 
NR --0 . 01022*0 . 5*RHO*L**4 
NV =-0. 00769*0. 5*RHO*L**3 
NDRS =-0.337*0. 02345*0. 5*RHO*L**3 
NDRB = + 0.283*0.02345*0. 5*RHO*L** 3 

WRITE (*,1001) 

READ (*,*) TCMI N , TCMAX , I TC 

WRITE (*,1002) 

READ (*,*) XDMIN.XDMAX, IXD 

XDMIN=XDMIN*L 

XDMAX-XDMAX *L 

WRITE (*,1003) 

READ ( * , * ) U 

WRITE ( * , 1100) 

READ ( * , * ) TL 

DH =( IZ-NRDOT) * ( MASS-YVDOT) - 
& ( MASS *XG- YRDOT ) * ( MAS S * XG-NVDOT ) 

AA1 1 - ( ( IZ-NRDOT) * YV- ( MASS *XG-YRDOT) *NV)/DH 
AAl 2- ( ( IZ-NRDOT) *(-MASS+YR)- 
& ( MASS *XG- YRDOT ) * ( -MASS * XG+NR ) )/DH 

AA2 1 - ( (MASS-YVDOT) *NV-(MASS*XG-NVDOT) *YV)/DH 
AA22- ( (MASS-YVDOT) * (-MASS *XG+NR)- 
& ( MASS * XG-NVDOT ) * ( -MASS + YR ) ) /DH 

BB 1 1 — ( ( IZ-NRDOT ) * YDRS- ( MASS * XG- YRDOT ) *NDRS ) /DH 
B B 1 2— ( ( IZ-NRDOT) * YDRB- (MASS* XG-YRDOT) *NDRB)/DH 
BB 2 1 - ( ( MASS-YVDOT ) * NDRS- ( MASS * XG-NVDOT ) * YDRS )/DH 
BB22= ( ( MASS-YVDOT) * NDRB- ( MASS *XG-NVDOT) * YDRB )/DH 

WRITE (*,1004) 
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c 



READ (*,*) 



RATIO 



C 



C 



C 



C 



C 



C 



BBl=BBll+RATIO*BBl2 
BB2=BB21+ RA TIO*BB22 

EPS =1 . D-5 
ILMAX=1500 

DO 1 1=1, ITC 

WRITE (*,2001) I, ITC 

TC=TCMIN+ ( 1-1 ) * ( TCMAX-TCMIN )/( ITC-1 ) 

TCD=TC* L/U 
POLEC-1 . 0/TCD 
ADl = 3 . 0*POLEC 
AD2=3 . 0*POLEC**2 
AD3=POLEC**3 
Al=BBl *U*U 
Bl=BB2*U*U 

Cl = - ADl - ( AAl 1 +AA2 2 ) *U 

A2= ( BB1 * AA2 2-BB2 *AAl 2 ) *U** 3 

B2 = ( BB2 * AAl 1 -BBl *AA2 1 ) *U* * 3 

Kl=AD3/( ( BB 2 *AA1 1-BB 1 *AA2 1 ) *U * * 3 ) 

C2=AD2-( AAl 1 *AA2 2-AA1 2 *AA2 1 )*U**2 + BB2*U*U*Kl 
K2=( Cl *B2-C2*Bl )/(Al*B2-A2*Bl) 

K3=(C2*Al-Cl*A2)/(Al *B2-A2 *B1 ) 

DO 2 J= 1 , IXD 

XD=XDMIN+( J-l ) * ( XDMAX-XDMIN )/( IXD-1 ) 

CALL L INEAR ( Kl , K2 , K 3 , U , XD , AAl 1 , AAl 2 , AA21 , AA22 , BB 1 , BB 2 , A , TL ) 
CALL RG(4,4,A,WR,WI,0,ZZZ,IV1, FV1 , I ERR) 

CALL DSTABL( DEOS ,WR, WI , FREQ) 

IF ( J .GT. 1 ) GO TO 10 

DEOSOO=DEOS 

XDOO =XD 

LL=0 

GO TO 2 

10 DEOSNN=DEOS 

XDNN =XD 
PR = DEOSNN *DEOSOO 
IF ( PR.GT. 0 . DO ) GO TO 3 
LL=LL+1 

IF (LL.GT.3) STOP 1000 
I L = 0 

XDO=XDOO 
XDN=XDNN 
DEOSO=DEOSOO 
DEOS N=D EOS NN 
6 XDL=XDO 

XDR=XDN 
DEOSL=DEOSO 
DEOS R=DEOSN 
XD= ( XDL+XDR )/2 . DO 

CALL LIN EAR ( Kl , K 2 , K 3 , U , XD , AAl 1 , AAl 2, AA21 ,AA22 , BBl,BB2,A,TL) 
CALL RG (4,4,A,WR,WI,0,ZZZ,IV1, FV1 , I ERR ) 

CALL DSTABLf DEOS ,WR,WI , FREQ) 

DEOSM=DEOS 

XDM=XD 
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PRL=DEOSL*DEOSM 

PRR-DEOSR'-DEOSM 

IF ( PRL .GT. 0 . DO ) GO TO 5 

XDO=XDL 

XDN=XDM 

DEOSO=DEOSL 

DEOSN-DEOSM 

I L= I L+l 

IF ( IL.GT. ILMAX) STOP 3100 
DI F=DABS ( XDL-XDM ) 

IF (DIF. GT . EPS ) GO TO 6 

XD=XDM 

GO TO 4 

5 IF (PRR.GT.0.D0) STOP 3200 

XDO=XDM 
XDN=XDR 
DEOSO=DEOSM 
DEOSN-DEOSR 
IL-IL+1 

IF ( IL.GT. ILMAX) STOP 3100 
D I F=DABS ( XDM-XDR ) 

IF (DIF. GT . EPS ) GO TO 6 
XD=XDM 

4 LLL = 1 0 + LL 

WRITE ( LLL , * ) XD/L , TC 
3 XDOO=XDNN 

DEOSOO=DEOSNN 
2 CONTINUE 
1 CONTINUE 



1001 


FORMAT 


( - 


ENTER 


MIN, 


MAX, 


AND 


INCREMENTS 


OF 


TC' ) 


1002 


FORMAT 


( ' 


ENTER 


MIN, 


MAX, 


AND 


INCREMENTS 


OF 


XD' ) 


1003 


FORMAT 


( ' 


ENTER 


U' ) 












1004 


FORMAT 


( ' 


ENTER 


BOW/STERN 


RUDDER RATIO' ) 






1100 


FORMAT 


(' 


ENTER 


TIME 


LAG 


TL' ) 








2001 


FORMAT 


(215) 















END 

C 

SUBROUTINE DSTABL ( DBOS , WR , WI , OMEGA) 

IMPLICIT DOUBLE PRECISION (A-H,0-Z) 

DIMENSION WR( 4 ) , WI ( 4 ) 

DEOS=-l . OD+20 
DO 1 1=1,4 

IF (WR( I ) .LT.DEOS) GO TO 1 
DEOS=WR ( I ) 

I J = I 

1 CONTINUE 

OMEGA=WI ( I J ) 

OMEGA=DABS ( OMEGA ) 

RETURN 

END 

C 

SUBROUTINE L INEAR( K1 , K2 , K3 , U , XD , AA1 1 , AA1 2 , AA2 1 , AA2 2 , BBl , BB2 , A , TL ) 
IMPLICIT DOUBLE PRECISION (A-H.O-Z) 

DOUBLE PRECISION K1,K2,K3 
DIMENSION A ( 4 , 4 ) 

A( 1 , 1 ) = 0 . 0D0 
A ( 1 , 2 ) =0 . 0D0 
A( 1 , 3 ) = 1 . 0D0 
A ( 1 , 4 ) = 0 . 0D0 
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A ( 2 , 1 )“(BBl*U*U*Kl-BBl*U*U*U*Kl *TL/XD )/ 

& ( 1-BBl *U*U*K1 *TL*TL/( 2 ,DO*XD) ) 

A ( 2 , 2 ) = ( AAl l*U+BBl*U*U*K2-BBl*U*U*Kl *TL/XD ) / 

& ( 1-BBl *U*U*Kl*TL*TL/( 2 . DO*XD) ) 

A ( 2 , 3 )=( AAl 2*U+BB1*U*U*K3+BB1 * U *U*U* Kl * TL* TL/( 2 . DO*XD) )/ 
& ( l-BBl*U*U*Kl *TL*TL/( 2 . DO*XD) ) 

A( 2, 4 )=( BBl*U*U*Kl/XD)/ 

& ( l-BBl*U*U*Kl * TL* TL/( 2 . DO * XD ) ) 

A ( 3 , 1 )■*( BB2*U*U*K1-BB2*U*U*U*K1 *TL/XD) + 

& ( BB2*U*U*Kl*TL*TL/( 2 . D0*XD) ) * 

& ( BBl*U*U*Kl-BBl*U*U*U*Kl *TL/XD ) / 

& ( l-BBl*U*U*Kl*TL*TL/( 2 ,D0*XD) ) 

A ( 3 , 2) = (AA21*U + BB2*U*U*K2-BB2 *U*U* K 1 *TL/XD ) + 

& (BB2*U*U*Kl *TL*TL/( 2 . D0*XD) ) * 

& ( AAl l*U+BBl*U*U*K2-BBl*U*U*Kl * TL/XD ) / 

& ( l-BBl*U*U*Kl *TL*TL/( 2 . D0*XD) ) 

A( 3, 3 )=( AA22*U+BB2*U*U*K3+BB2*U*U*U*Kl*TL*TL/( 2.D0*XD) )+ 
& (BB2*U*U*Kl* TL *TL/ ( 2 . DO * XD ) ) * 

& ( AA12*U+BB1 *U*U*K3+BBl *U*U*U*K1 *TL*TL/( 2 . DO*XD) )/ 

& ( l-BBl*U*U*Kl *TL*TL/( 2 .D0*XD) ) 

A ( 3 , 4 ) = ( BB2*U*U*Kl/XD) + ( BB2*U*U*Kl *TL*TL/( 2 . D0*XD) ) * 

& ( BB1 *U*U* Kl/XD ) / 

& ( l-BBl*U*U*Kl *TL*TL/( 2 . DO*XD) ) 

A ( 4 / 1 ) =U 

A ( 4 , 2 ) « 1 . ODO 

A ( 4 , 3 )=0. ODO 

A ( 4 , 4 ) “0 . ODO 

RETURN 

END 
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APPENDIX C 



PROGRAM THES I S 3 . FOR (Time Delay-3rd Order Approx) 
BIFURCATION ANALYSIS 

IMPLICIT DOUBLE PRECISION (A-H.O-Z) 

DOUBLE PRECISION K 1 , K2 , K3 , L , NR , NV , NDRS , NDRB , I Z , MASS , 

& NRDOT , NVDOT 

DIMENSION A( 5 , 5 ) , B ( 5 , 5 ) , FVl ( 5 ) , IVl ( 5 ) , ZZZ ( 5 , 5 ) , ALFR ( 5 ) , 
& ALF I ( 5 ) , BETA (5),WR(5),WI(5) 

OPEN (11, FILE-' BI FI .RES' , STATUS- ' NEW ' ) 

OPEN ( 12, FILE- 'BIF2. RES' , STATUS- ' NEW ' ) 

OPEN <13,FILE='BIF3.RES' , STATUS- ' NEW' ) 

Vehicle Parameters: 

WEIGHT-435.0 
IZ =45.0 
L =7.3 

RHO -1.94 
G -32.2 

XG -0.0104 
MASS -WEIGHT/G 

YRDOT — 0. 00178*0. 5*RHO*L**4 
YVDOT =-0 . 03430*0 . 5*RHO*L**3 
YR =+0 . 01187*0 . 5*RHO*L* *3 
YV --0 . 03896*0 . 5*RHO*L**2 
YDRS -+0 . 02345*0 . 5*RHO*L**2 
YDRB -+0 . 02345*0 . 5*RHO*L**2 
NRDOT --0. 00047*0. 5*RHO*L**5 
NVDOT --0. 00178*0. 5*RHO*L**4 
NR «-0 .01022*0 . 5*RHO*L**4 
NV --0 . 00769*0 . 5*RHO*L**3 
NDRS =-0.337*0 . 02 345 * 0 . 5 *RHO*L* * 3 
NDRB = + 0.283*0. 02345*0 . 5*RHO*L** 3 

WRITE (*,1001) 

READ (*,*) TCMIN , TCMAX , I TC 

WRITE (*,1002) 

READ (*,*) XDMIN , XDMAX , I XD 

XDMIN-XDMIN* L 
XDMAX- XDMAX * L 
WRITE (*,1003) 

READ ( * , * ) U 
WRITE( *,1100) 

READ ( * , * ) TL 

DH = ( I Z -NRDOT ) * ( MASS- YVDOT ) - 
& ( MASS * XG- YRDOT ) * ( MASS *XG-NVDOT ) 

AA1 1 - ( ( I Z-NRDOT ) * YV- ( MASS *XG- YRDOT ) *NV )/DH 
AA1 2- ( ( I Z-NRDOT) * ( -MASS+YR)- 
& ( MASS *XG- YRDOT ) * ( -MASS *XG+NR ) )/DH 

AA2 1 = ( ( MASS-YVDOT ) *NV- ( MASS *XG-NVDOT ) * YV )/DH 
AA22- ( (MASS-YVDOT) * ( -MASS *XG+NR ) - 
& ( MASS *XG-NVDOT ) * ( -MASS + YR ) )/DH 

BB 1 1 — ( ( I Z-NRDOT) * YDRS- ( MASS *XG- YRDOT ) *NDRS )/DH 
BBl 2- ( ( I Z-NRDOT ) * YDRB- (MASS* XG-YRDOT) *NDRB)/DH 
BB2 1 — ( (MASS-YVDOT ) *NDRS- ( MASS * XG-NVDOT ) * YDRS ) /DH 
BB22- ( (MASS-YVDOT) *NDRB- ( MASS * XG-NVDOT ) * YDRB ) /DH 
C 

WRITE (*,1004) 
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c 



READ (*,*) 



RATIO 



C 



C 



C 



C 



C 



BBl=BBll+ RATIO * BB 1 2 
BB2-BB21+RATIO*BB22 

EPS “l.D-5 
I LMAX= 1500 

DO 1 1=1, ITC 

WRITE (*,2001) I, ITC 

TC=TCMIN+( 1-1 ) * ( TCMAX-TCM I N ) / ( ITC-1 ) 

TCD=TC * L/U 
POLEC-1 . 0/TCD 
ADI-3 . 0*POLEC 
AD2-3 ,0*POLEC**2 
AD3=POLEC**3 
A1=BB1 *U*U 
Bl-BB2*U*U 

Cl=-ADl-( AA11+AA22 ) *U 

A2= ( BB1 ‘AA22-BB2 * AA1 2 ) *U**3 

B 2 = ( BB2 * AAl 1 -BB 1 * AA2 1 )*0**3 

Kl = AD3/( ( BB2*AAll-BBl *AA21 ) *U* *3 ) 

C2=AD2- ( AAl 1 * AA2 2 -AAl 2 * AA2 1 )*U**2 + BB2*U*U*K1 
K2 = (Cl*B2-C2*Bl )/( Al * B2-A2 * B1 ) 

K3=(c2*Al-Cl*A2 )/( Al *B2-A2*B1 ) 

DO 2 J = 1 , IXD 

XD=XDMIN+ ( J-l ) * ( XDMAX-XDMIN )/ ( IXD-1 ) 

CALL LINEAR ( Kl , K 2 , K 3 , U , XD , AAl 1 , AAl 2 , AA2 1 , AA22 , BB 1 , BB2 , A , B , TL ) 
CALL RGG ( 5, 5, A, B , ALFR, ALFI , BETA, 0 , ZZZ, I ERR ) 

DO 11 IJE=1,5 

WR ( I J E ) =ALFR ( 1 0 E ) /BETA ( IJE) 

WI ( IJE) =ALF I ( IJE )/BETA ( IJE) 

11 CONTINUE 

CALL DSTABL ( DEOS,WR,WI , FREQ) 

IF ( J .GT. 1 ) GO TO 10 

DEOSOO=DEOS 

XDOO =XD 

LL = 0 

GO TO 2 

10 DEOSNN=DEOS 

XDNN =XD 
PR=DEOSNN*DEOSOO 
IF ( PR.GT. 0 . DO ) GO TO 3 
LL-LL+1 

IF (LL.GT.3) STOP 1000 
I L = 0 

XDO=XDOO 
XDN=XDNN 
DEOSO=DEOSOO 
DEOSN-DEOSNN 
6 XDL-XDO 

XDR=XDN 
DEOSL-DEOSO 
DEOSR=DEOSN 
XD=(XDL+XDR)/2 . DO 

CALL LI NEAR ( Kl , K2 , K 3 , U , XD , AAl 1 , AAl 2 , AA2 1 , AA2 2 , BB 1 , BB 2 , A , B , TL ) 
CALL RGG ( 5, 5, A, B, ALFR, ALFI , BETA, 0, ZZZ, I ERR) 
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12 

C 

5 

4 

3 

2 

1 

C 

1001 

1002 

1003 

1004 

1100 

2001 

C 

1 

c 



DO 12 IJE-1,5 

WR ( IJ E ) ~ALFR ( I J E )/BETA( I JE) 
WI( IJE)«ALFI( I J E ) /BETA ( IJE) 
CONTINUE 

CALL DSTABL ( DEOS , WR , WI , FREQ ) 

DEOSM-DEOS 

XDM-XD 

PRL-DEOSL * DEOSM 

PRR-DEOSR * DEOS M 

IF (PRL.GT.0.D0) GO TO 5 

XDO-XDL 

XDN-XDM 

DEOSO-DEOSL 

DEOSN-DEOSM 

IL-IL+1 

IF ( IL.GT. ILMAX) STOP 3100 
DIF=DABS( XDL-XDM) 

IF (DIF. GT . EPS ) GO TO 6 

XD-XDM 

GO TO 4 

IF ( PRR.GT.0.D0) STOP 3200 

XDO=XDM 

XDN=XDR 

DEOSO-DEOSM 

DEOSN-DEOSR 

I L= I L+ 1 

IF ( IL.GT. ILMAX) STOP 3100 
DI F=DABS ( XDM-XDR ) 

IF ( DIF.GT.EPS) GO TO 6 

XD=XDM 

LLL=* 1 0 + LL 

WRITE ( LLL , * ) XD/L , TC 
XDOO-XDNN 
DEOSOO-DEOSNN 
CONTINUE 
CONTINUE 



FORMAT 


( ' 


ENTER 


MIN, 


MAX, 


AND 


INCREMENTS 


OF 


TC' ) 


FORMAT 


( ' 


ENTER 


MIN, 


MAX, 


AND 


INCREMENTS 


OF 


XD' ) 


FORMAT 


( ' 


ENTER 


U' ) 












FORMAT 


(' 


ENTER 


BOW/STERN 


RUDDER RATIO' ) 






FORMAT 


( ' 


ENTER 


TIME 


LAG 


TL' ) 









FORMAT (215) 

END 

SUBROUTINE DSTABL ( DEOS,WR,WI , OMEGA) 
IMPLICIT DOUBLE PRECISION (A-H,0-Z) 
DIMENSION WR ( 5 ) , WI ( 5 ) 

DEOS**- 1 . 0D+20 
DO 1 1-1,5 

IF (WR( I ) . LT.DEOS ) GO TO 1 
DEOS-WR ( I ) 

I J«I 

CONTINUE 
OMEGA-WI ( I J ) 

OMEGA=DABS( OMEGA) 

RETURN 

END 
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SUBROUTINE LINEAR! Kl , K2 , K 3 , U , XD , AAl 1 , AA12 , AA21 , AA22 , BBl , BB2 , 
& A , B , TL ) 

IMPLICIT DOUBLE PRECISION (A-H,0-Z) 

DOUBLE PRECISION Kl,K2,K3 
DIMENSION A( 5, 5 ) , B( 5, 5 ) 

A ( I,1)-0.0D0 
A ( 1,2) “0 . 0D0 
A( 1 , 3 ) -1 . 0D0 
A( 1 , 4 )=0 . 0D0 
A( 1 , 5 ) =0 . OdO 
A( 2 , 1 )-0 . ODO 
A( 2 , 2 ) =0 . ODO 
A(2, 3 ) =0 . ODO 
A ( 2 , 4 ) -0 . OdO 
A( 2 , 5 )=1 . ODO 

A( 3, 1 )-(BBl*U*U*Kl-BBl*U*U*U*Kl*TL/XD) 

A ( 3 , 2 )=( AAll *U+BBl *U*U*K2-BBl*U*U*Kl *TL/XD) 

A ( 3 , 3 )»( AA12*U+BB1*U*U*K3+BB1 * U*U * U* Kl * TL*TL/( 2 . DO*XD) ) 

A ( 3 , 4 ) =■( BBl*U*U*Kl/XD ) 

A ( 3 , 5 ) = ( - 1 . DO + BBl *U*U*Kl*TL*TL/( 2 .D0*XD) ) 

A ( 4 , 1 ) =U 
A ( 4 , 2 ) = 1 . ODO 
A( 4, 3)-0.0D0 
A( 4 , 4 )-0 . ODO 
A( 4 , 5 )=0 . OdO 

A( 5, 1 ) = (BB2*U*U*Kl-BB2*U*U*U*Kl*TL/XD) 

A( 5, 2 ) = ( AA2 l*U+BB2*U*U*K2-BB2*U*U*Kl *TL/XD) 

A( 5, 3 )-( AA22*U+BB2*U*U*K3+BB2*U*U*U*K1 *TL*TL/( 2 . D0*XD) ) 

A ( 5 , 4 )-( BB2*U*U*K1/XD) 

A( 5, 5 ) - ( BB2*U*U*Kl*TL*TL/( 2 . DO *XD ) ) 

b ( l , l ) - l . O dO 
B( 1 , 2 ) -0 . ODO 
B( 1 , 3 )=0 . ODO 

b ( 1 , 4 ) =0 . OdO 
b ( 1 , 5 ) = 0 . ODO 
B ( 2 , 1 ) =0 . ODO 
B( 2, 2)-l .ODO 
B( 2, 3 ) =0 . OdO 
B ( 2 , 4 ) “0 . ODO 
B( 2, 5)=0 . ODO 
B( 3 , 1 ) =0 . ODO 
B( 3, 2)=0.0D0 

B( 3 , 3 )=( BBl *U*U*U*Kl*TL*TL*TL/( 6 . D0*XD) ) 

B( 3, 4 )=0 . ODO 

B( 3 , 5 )-( BBl *U*U*Kl *TL*TL*TL/( 6 . D0*XD) ) 

B( 4 , 1 ) =0 . ODO 
B ( 4 , 2 ) =0 . OdO 
B( 4, 3 )-0 . ODO 
B( 4 , 4 )-l .ODO 
B( 4, 5)»0 .OdO 
B ( 5 , 1 ) “0 . ODO 
B ( 5 , 2 ) =0 . ODO 

B( 5, 3)-( 1 . 0D0+BB2*U*U*U*Kl*TL*TL*TL/( 6 . D0*XD ) ) 

B( 5, 4 ) -0 . ODO 

B( 5 , 5 )-( BB2*U*U*Kl*TL*TL*TL/( 6 . D0*XD) ) 

RETURN 

END 
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APPENDIX D 



PROGRAM THESISN.FOR (Time Delay - Newton's Iteration Method) 
IMPLICIT DOUBLE PRECISION (A-H.O-Z) 

DOUBLE PRECISION K1 , K2 , K 3 , L , NR , NV , NDRS , NDRB , I Z , MASS , 

& NRDOT , NVDOT 

DIMENSION A ( 4 , 4 ) , FV1 ( 4 ) , IVl ( 4 ) , ZZZ ( 4 , 4 ) , WR( 4 ) , WI ( 4 ) 

OPEN ( 11 , FILE-' BI Fl . RES' , STATUS- 'NEW' ) 

Vehicle Parameters: 

WEIGHT-435.0 
IZ -45.0 

L —7 . 3 

RHO =1.94 

G -32.2 

XG =0.0104 

MASS -WEIGHT/G 
C 

YRDOT =-0 . 00178*0 . 5*RHO*L**4 
YVDOT =-0 . 03430*0 . 5*RHO*L** 3 
YR =+0. 01187*0. 5*RHO*L**3 

YV =-0 . 03896*0 . 5*RHO*L* *2 

YDRS =+0.02345*0. 5*RHO*L**2 
YDRB - + 0. 02345*0. 5*RHO*L**2 
NRDOT =-0. 00047*0. 5*RHO*L**5 
NVDOT =-0 . 00178*0 . 5*RHO*L**4 
NR =-0 . 01022*0 . 5*RHO*L**4 

NV =-0. 00769*0. 5*RHO*L**3 

NDRS =-0. 337*0.02 345*0.5*RHO*L**3 
NDRB =+0.283*0.02345*0 . 5*RHO*L** 3 
C 

WRITE (*,1001) 

READ (*,*) TCMIN, TCMAX , ITC 

WRITE (*,1003) 

READ ( * , * ) U 

C 

DH = ( IZ-NRDOT ) * ( MASS-YVDOT ) - 

& ( MASS *XG- YRDOT ) * ( MASS*XG-NVDOT ) 

AA1 1 = ( ( IZ-NRDOT ) * YV- ( MASS* XG- YRDOT ) *NV ) /DH 
AA 1 2 = ( ( IZ-NRDOT) * ( -MASS+YR)- 
& ( MASS *XG- YRDOT ) * ( -MASS*XG+NR ) )/DH 

AA2 1 = ( (MASS-YVDOT) *NV- ( MASS *XG-NVDOT ) *YV)/DH 
AA22=( ( MASS-YVDOT) * ( -MASS *XG+NR ) - 
S ( MASS * XG- NVDOT ) * ( -MASS + YR ) )/DH 

BB 1 1 — ( ( IZ-NRDOT) * YDRS- ( MASS*XG- YRDOT ) *NDRS)/DH 
BB 1 2= ( ( IZ-NRDOT ) * YDRB- ( MAS S*XG- YRDOT ) *NDRB ) /DH 
BB21=( (MASS-YVDOT) *NDRS- ( MASS *XG-NVDOT ) * YDRS ) /DH 
BB22=( (MASS-YVDOT) *NDRB- ( MASS *XG-NVDOT ) * YDRB )/DH 
C 

WRITE (*,1004) 

READ (*,*) RATIO 

C 

BBl = BBll +RATIO* BB 1 2 
BB2=BB21+RATIO*BB22 
C 

WRITE (*,1005) 

READ (*,*) TL 

C 

EPS= 1 . D- 1 0 
ITMAX-1000 
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c 

DO 1 1=1 , ITC 

WRITE (*,2001) I, ITC 

TC-TCMI N+ ( 1-1 ) * (TCMAX-TCMIN)/( ITC-1 ) 

TCD=TC * L/U 
POLEC-1 . 0/TCD 
AD 1 = 3 . 0 * POL EC 
AD2=3 . 0*POLEC* *2 
AD3=POLEC**3 
Al=BBl*U*U 
B1 = BB2 *U*U 

Cl=-ADl-( AA1 1+AA22 ) *U 
A2= ( BBl * AA22-BB2 * AAl 2 ) *U**3 
B2= ( BB2 * AAl 1-BB 1 * AA2 1 )*U**3 
K 1 =AD 3/ ( ( BB2 * AA 1 1 — BB 1 * AA2 1 ) *U* * 3 ) 

C2=AD2-( AAl 1 * AA22-AA1 2 * AA2 1 ) *U**2 + BB2*U*U*K1 
K2=(C1*B2-C2*B1 )/< A1*B2-A2*B1 ) 

K3=(C2*Al-Cl*A2)/(Al *B2-A2 * B 1 ) 

Al- ( AAl 1 * BB2-AA2 1*BB1 ) *U*U*U*Kl 

A2= ( AAl 1 * AA2 2-AAl 2 * AA2 1 ) *U*U+ ( AAl 1 *BB2-AA21 *BBl )*U*U*U*K3+ 

& (AA22*BB1-AA12*BB2) *U*U*U*K2-BB2*U*U*K1 

A3=- ( AA11+AA22 )*U-(BBl*K2+BB2*K3)*U*U 
B0= (AAll*BB2-AA21*BBl ) *U*U*U*U*K1 
B 1 = - ( AAl 2 * BB2-AA2 2 * BBl )*U*U*U*Kl-BB2*U*U*U*Kl 
B2 — BBl*U*U*Kl 

COMPUTE INITIAL APPROXIMATION FOR OMEGA (TL-0) 

AQ=B2 * A 3-B 1 
BQ=Bl *A2-B2 *A1-B0 *A3 
CQ=B0*Al 

DET=BQ* *2-4 .D0*AQ*CQ 
IF ( DET. LT. 0 . DO ) GO TO 1 
ZT1 = ( -BQ+DSQRT ( DET ) )/(2.0D0*AQ) 

ZT2= ( -BQ-DSQRT ( DET ) )/(2.0D0*AQ) 

IF ( ZTl .GE . 0 . DO ) OM=DSQRT( ZTl ) 

IF ( ZT2 .GE . 0 .DO ) OM=DSQRT( ZT2 ) 

ALPHAl=AQ 
ALPHA2=BQ 
ALPHA3=CQ 
B E TA 1 = - B 2 

BETA2 = B0 + B2 *A2-B1*A3 
BETA3=B1 *A1-B0*A2 

COMPUTE EXACT OMEGA USING NEWTON'S ITERATION 
OMOLD=OM 

3 F=(BETAl*OMOLD**5+BETA2*OMOLD**3+BETA3*OMOLD) *DSIN(OMOLD*TL) 
& + ( ALPHAl *OMOLD* * 4 +ALPHA2 *OMOLD* * 2+ALPHA3 ) * DCOS ( OMOLD * TL ) 

FPRIMEM 5 . DO * BETA 1 * OMOLD** 4 + 3 . DO * BETA 2 * OMOLD* * 2 + BET A3 ) 

& *DSIN(OMOLD*TL)+(BETAl*OMOLD**5+BETA2*OMOLD**3+BETA3*OMOLD) 
& *TL*DCOS (OMOLD * TL)+ ( 4 . DO * ALPHAl * OMOLD * * 3 + 2 . DO *ALPHA2 * OMOLD ) 

& * DCOS ( OMOLD *TL)-(ALPHAl * OMOLD* * 4+ALPHA2 * OMOLD* * 2 + ALPHA 3 ) 

& *TL*DSIN( OMOLD* TL ) 

IF (FPRIME.EQ.0.D0) STOP 1112 
IF (F.EQ.0.D0) GO TO 2 
OMNEW=OMOLD-F/FPRI ME 
OMDI F=DABS ( OMNEW-OMOLD ) 

IF (OMDIF.LE.EPS) GO TO 2 
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OMOLD-OMNEW 

IT-IT+1 

IF < IT.GT. ITMAX) STOP 1111 
GO TO 3 
C 

2 OM-OMNEW 

XDNUM-DSQRT ( (BO-B2*OM**2)**2+(B1*OM)**2) 
XDDEN=OM*DSQRT< ( Al -A3 *OM * * 2 ) * * 2+ ( A2 *OM-OM* * 3 ) * * 2 ) 
XD=XDNUM/XDDEN 
XD=XD/L 
C 

WRITE (11,*) XD,TC,OM 
C 

1 CONTINUE 



1001 FORMAT 

1003 FORMAT 

1004 FORMAT 

1005 FORMAT 
2001 FORMAT 

END 



(' ENTER MIN, MAX, AND INCREMENTS OF TC ' ) 
< ' ENTER U' ) 

( ' ENTER BOW/STERN RUDDER RATIO' ) 

< ' ENTER TIME LAG' ) 

( 215 ) 
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